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Abstract The stellar evolution code YREC is outlined with 
emphasis on its applications to helio- and asteroseismology. 
The procedure for calculating calibrated solar and stellar 
models is described. Other features of the code such as a 
non-local treatment of convective core overshoot, and the 
implementation of a parametrized description of turbulence 
in stellar models, are considered in some detail. The code 
has been extensively used for other astrophysical applica- 
tions, some of which are briefly mentioned at the end of the 
paper. 

Keywords methoods: numerical • stars: evolution • stars: 
interior • convection 
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O , 1 Introduction 



The aim of this paper is to provide an overview of the Yale 
Rotating Stellar Evolution Code (YREC), as it has been ap- 
plied in the last few years to research in helio- and astero- 
seismology. Although YREC contains extensions to model 
the effects of rotation in an oblate coordinate system, we de- 
scribe here the "non-rotating" version. 
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In addition to a general description, we shall emphasize 
three features of the code which have been implemented be- 
cause of their special relevance to seismology. The first fea- 
ture is the procedure utilized for the automatic calculation of 
calibrated solar and stellar models whose pulsational prop- 
erties are to be investigated. The second feature is the treat- 
ment of convective core overshoot. Finally, the third feature 
is the implementation in stellar models of the effects of tur- 
bulence on the structure of the surface layers of stars with 
a convective envelope. The parametrization of turbulence to 
one dimension is based on three-dimensional radiative hy- 
drodynamical (3D HRD) simulations of the highly supera- 
diabatic layer (SAL) in the atmosphere. The interaction of 
turbulent convection and radiation in these thin transition 
regions is poorly known. Oscillation frequencies are sensi- 
tively affected by the structure of transition regions between 
radiative and convective layers. Seismology thus offers a 
unique opportunity to explore a long standing problem in 
stellar physics. 

Like most stellar evolution codes, YREC is a continu- 
ously evolving research tool to which many have contributed. 
As a result, different versions of YREC are in use at sev- 
eral institutions, which have been applied to a variety of re- 
search purposes. Some of the most significant applications 
of YREC are listed in the text and at the end of this paper 
(see Sect. |T3T>. The rotating vers ion of YREC, originally de- 
veloped by iPinsonneaultl (Il988l) . in cludes a 1.5D treatment 
of rota tion, extending the work of [Kippenhahn & Thomas! 
(119701) and lEndal & Sofial dl98l[) . and using the formalism 
of lLawl (1 19801) . A 2D version of YREC has also recently 
been implemented, specifically to ad dress some fun damen- 
tal aspects of solar magnetic activity (ILi et al.ll200q) . 

Sect. [2] outlines the numerical scheme adopted to solve 
the classical differential equations of stellar structure and 
evolution. The treatment of the boundary conditions, of spe- 
cial importance for seismology, are described in Sect.[3j The 
constitutive physics, i.e. the equation of state and radiative 
and conductive opacities, are reviewed in Sect. |H and the 
nuclear processes are described in Sect. [6] Stellar physics 
topics such as superadiabatic convection, element diffusion, 
convective core overshoot, and turbulence in the outer lay- 
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ers, all of which also have important seismological signa- 
tures, are covered in Sect. [321 Sect. [5] Sect. [9] and Sect.flOl 
respectively. The operation of the code is described in Sect.[7] 
with emphasis on helio- and asteroseismic applications. Seis- 
mic diagnostics applications are described in Sect. [TTJ The 
role played by YREC in the research on solar neutrinos and 
helio-seismology is summarized in Sect. [T2j Studies of ad- 
vanced evolutionary phases and applications to stellar popu- 
lation studies are listed in Sect. [13] 



2 Henyey code 

The four first-order simultaneous equations of stellar struc- 
ture are well- known, and have bee n frequently discussed in 
the literature dSchwarzschildlll958h . YREC uses mass as the 
independent variable in the formulation of the equations (La- 
grangian formulation). The problem is a two-point boundary 
value problem, with boundary conditions at the center and 
at the surface of the model. A relaxation technique, based 
on a finite difference approximation, is used. The method, 



of hydrogen it begins burning hydrogen in a shell. Initially 
the shell is almost 0.2M Q thick for a 1M© star but the shell 
quickly narrows to only O.OOIM^ as the star evolves up the 
giant branch, thinning to O.OOOIM© at the point of helium 
flash. Because of the high temperature dependence of he- 
lium burning, helium burning shells are even thinner than 
hydrogen burning shells. YREC will add or remove shells 
according to the size of gradients in structure (i.e., pressure, 
temperature, and composition) and gradients in luminosity, 
as well as the size of Henyey corrections applied during the 
iteration procedure. The code keeps track of physically real 
discontinuities so that they are not smoothed during the re- 
zoning process. Interpolation is linear. Our own testing has 
shown that using higher order methods such as oscillatory 
spline interpolation introduces numerical oscillations near 
the tip of the giant branch. 



2.2 Time steps 



first ap plied to the stellar structure problem bv lHenvev et all Th e models are advanced in time through two terms in the 

energy equation, the nuclear energy term (Sect. [6J|, and the 
arson & DemarffiHl rate °f cnan g e of entropy due to contraction or expan 



first ap p 

dl959h . is known as the Henyey method. Usef ul descriptions 



of the Henyey method are given in the paper by 



dl964l) and the book bv lKippenhahn & Weigertldl990l) . Spe 



cific details about the numerical procedures in the YREC 
implem entation can be found in the Appendices of IPratherl 
dl976h . which d escribe an earl ier Henyey code on which 
YREC is based. I Prattler! d 19761) also provides information 
about the treatment of the constitutive physics, although most 
of the physics details have been updated since then. 

In the Henyey method, the model star is divided into n 
concentric shells by means of n + 1 suitably chosen values 
of the independent variable (mass), or points, in the interval 
defined by the innermost point (near the center) and the out- 
ermost point, which is specified by the user and located at 
the base of the envelope integrations. The four differential 
equations are replaced in each shell by approximating dif- 
ference equations relating the values of the dependent vari- 
ables at adjacent points. There are four dependent variables 
in each of the n shells, providing a set of 4(n + 1) linear 
equations which, together with two boundary conditions at 
the center and two at the surface, can be solved to determine 
approximate corrections to the 4(n + 1 ) dependent variables, 
starting from a first approximation model. The set of simul- 
taneous equation is solved by iteration until the corrections 
in each variable satisfy a specified convergence limit. 



sion during evolution. Special care is taken to prese rve nu- 
merical accuracy for small time steps dPratherl dl976t) l 

One can either specify the time step or have YREC au- 
tomatically determine the optimum time step during evolu- 
tion. When producing accurately calibrated solar models, 
to maintain numerical consistency it helps to specify the 
time step interval. In most other situations it is best to let 
YREC determine the time step based on user specified con- 
vergence tolerance criteria. During nuclear burning phases 
of evolution YREC will guess the time step based on the 
rate at which hydrogen and helium (if applicable), are being 
consumed in each shell of the model. During gravitational 
contraction phases of evolution YREC will control the time 
steps by monitoring the change in temperature, pressure, and 
luminosity from one model to the next. During helium flash 
if the model fails to converge during a evolutionary time 
step, YREC is also able reduce the time step by a user spec- 
ified factor and redo the evolutionary step. More details re- 
garding the operation of YREC can be found in Sect. [7] 



3 Boundary conditions 



2.1 Shell redistribution 

The shells are distributed so as to optimize numerical accu- 
racy and efficiency. In order to follow the evolution from the 
earliest gravitational contraction phase all the way to the hy- 
drogen and helium shell burning phases, it is necessary to 
redistribute the shells in the model. This is especially criti- 
cal during shell burning. After a star exhausts its core supply 



3.1 Center 

The two inner boundary conditions constrain the values of 
the radial distance and luminosity variables at the innermost 
mass shell. Because of the false singularity at the center, the 
innermost point is not at the very center, but in a shell chosen 
close to the center. Note that in order to preserve accuracy, 
special care must be taken with the position of the innermost 
shell, especially in pulsation calculations (see Sect. 17.31 ). 
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3.2 Surface 



4 Equation of state and opacities 



The outer boundary conditions depend on the structure near 
the surface. Because a model of the outer layers depends on 
the global properties of the star, i.e. its surface gravity log g 
and effective temperature r e ff, the problem is implicit. In or- 
der to specify the surface boundary conditions, which relate 
the variations of the pressure and temperature variables to 
the total luminosity and radius of the star, three inward en- 
velope integrations are constructed. These envelope integra- 
tions are chosen so as to form a triangle in the theoretical 
HR-diagram (i.e. the log L/Lq vs. log T e g plane). 

The inward envelope integrations consist of two main 
parts. The outermost layers, starting at optical depth near 
T= 10" 10 , which are effectively isothermal at the start, are 
described by a gray radiative atmosphere specified by a T(z) 
relation and integrated to the appropriate value of T at which 
the temperature reaches T e ff (e.g. T = 2 /3 for the Edding t on 
approximation, t = 0.312156330 for fhe lKrishna Swamvl dl96 
atmosphere). This surface in the star is usually defined as the 
photosphere. 



As an alternative to the atmosphere integrations, more 
complex atmospheres from pre-computed libraries can be 
also used, such as those from lKuruczl dl998l) . 

Below the photosphere, all variables but the luminosity 
variable (which is held to be constant in the outer envelope) 
are integrated to a chosen value of the mass. The integra- 
tion is carried out using log P as the independent variable, 
to the value of the mass variable at which the surface bound- 
ary conditions for the interior are computed (the base of the 
envelope). The region which extends from this value of the 
mass to the innermost shell of the star constitutes the in- 
terior of the stellar model. In convectively unstable layers 
of the envelope (as determined by the local Schwarzschild 
criterion), the tempe rature gradient i s evalu ated according 
to the formalism of IStothers & Chin] d 19951) . which is de- 
signed to describe superadiabatic convection. It is in this 
region that the peak of the highly superadiabatic transition 
layer (SAL) is norma lly located (as it is in t he Sun). The 
main advantage of the IS tothers & Chinl dl995l) formalism is 
that by a suitable choice of parameters, it can be made to 
reproduce either the standard mixing length theory (MLT) 
Bohm- Vitensd 1 1 95 8l) or the theory of ICanuto & Mazzitellil 



YREC has been updated regularly so as to incorporate the 
latest research developments regarding equation of state and 
opacity in the stellar interior, while maintaining backward 
compatibility with earlier versions of the same. The current 
versio n uses the latest OPAL op acities ( Iglesias & Rogers! 
19961) and OPAL equation of state dRogers & Navfonovfcoda) . 
At low temperatures (log T < 4.1) opacities are obtained 
from the tables provided by lFerguson et al.l d2005l ). 

At each mass shell the EOS is obtained by interpolation 
from the standard tables. Since the EOS is weakly dependent 
on Z, we use only one set of tables at a fixed Z, obtained 
by the Z-interpolation routine provided with the OPAL EOS 
package. For models with metal diffusion, the value of Z at 
which the EOS is interpolated is chosen at a suitable inter- 
mediate value. The EOS quantities at the desired X, T and p 
are obtained by quadratic interpolation from the tables. The 
results and the derivatives are smoothed by mixing overlap- 
ing quadratics. For opacity, a four-point Lagrangian inter- 
olation scheme is used over a 4-dimensional grid of Z, X, 
T, and p. 



5 Diffusion 

The diffusion of chemical elements by gravitational settling 
and the rmal diffusion is im plemented following the prescrip- 
tion of iThoul et al.l d 19941) . Options in the code include no 
diffusion, helium diffusion onl y, or both Y and Z diffusion. 
The analytical fits provided by IThoul et all dl994h can also 
be used instead of the tabulated diffusion coefficients, to 
speed up the computations. 



6 Nuclear Reactions 

The nuclear reaction rates in conjunction with the corre- 
sponding energy release (Q-values) are important for the 
evolution of chemical species, the energy input from nuclear 
fusion reactions and for the neutrino fluxes. The reactions 
explicitly calculated in YREC are the following: 



19921) . sometimes called FST. In order to preserve continu- 



ity in the convective temperature grad i ent at the envelope- 
interior interface, the IStothers & Chir] dl995l) formalism is 
used to calculate the convective gradient both in the enve- 
lope and in the interior whenever superadiabaticity exceeds 
a preset value. 

Another feature of the envelope integration, described in 
more detail in Sect.QIJl includes a ID parametrization of the 
effects of turbulent pressure and turbulent kinetic energy in 
the outer layers. 
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The first five Equations (fl]|5]l contain the three alternative pp 
branches (ppl,pp2,pp3) all of which start with 3 He. Equa- 
tions ©-(O represent the primary CNO cycle, 1 Ob the sec- 
ondary cycle. The reaction of helium burning is given in 
Equation dill ), followed by the dominant a capture reac- 
tions d 1 2b - d 14b - The last two reactions are only important 
for the neutrino problem and can be neglected for the energy 
generation. As is implicitly shown in the nuclear reactions 
dTlll6l). all P -decay reactions are treated in the instantaneous 
approximation. In additi on, four branching ratios are defined 
dBahcall & Ulrichll 19881) : the fraction of 7 Be that is burned 
by electron capture ©, the fraction of 7 Be that is burned 
by proton capture (0), the fraction of 14 N that is burned 
via 14 N(p,a) 12 C and the fraction of 15 N that is burned via 

15 N(p,y) 16 0. 

The energy generation is calculated b y multiplying the 
rates b y the Q-values which are taken from lBahcall & Ulrichl 
dl988l Table 21). The standard reac tion rates imple mented 
are identical to the rates published in lBahcalll d 19891 ). 



6.1 NACRE 

YRE C provides the opti on of using the NACRB3 reaction 
rates dAngulo et al.ll999l) . In its present version, the Q-values 
from each reaction are not changed and are thus not identi- 
cal to the values published on the NACRE databasfl All 
relevant reaction rates that are provided by NACRE are in- 
cluded, i.e. those corresponding to Equations dTT3l I51ITTV 

The NACRE library lists the rate data in tabulated form 
and also provides fit-formulas, the latter of which are im- 
plemented. The fit-formulas are accurate by 3% - 25% com- 
pared to the tabulated data, with typical deviations of 10% - 
15%. 

Our tests for a standard solar model have shown that a 
calibrated standard model is not affected by the NACRE re- 
action rates. The largest differences are found in the neutrino 
flux of 8 B which differs by about 9%. This diff erence is com- 
fortab ly within other theoretical uncertainties dBahcall et alj 
2004). 



6.2 Light elements 

A switch permits keeping track of nuclear burning of the 
light elements 2 H, 6 Li, 7 Li and 9 Be at t he base of the cq n- 
vection zone in models of sun-like stars dDeliy annisll 1 990h . 



Nuclear Astrophysics Compilation of REaction Rates 
2 http://pntpm.ulb.ac.be/Nacre 



6.3 Neutrino losses 

Neutri no loss rates are taken from the monograph by lBahcaill 
dl989h . updated by subsequent private communications from 
the author. For advanced stages of stellar evolution, the neu- 
trino r ates from photo, pair and plasma sources from lltoh et alj 
d 19891 ) are included. 



7 Running YREC 

YREC can automatically calculate calibrated solar and stel- 
lar models. The user provides a complete set of constraints 
along with allowable parameter variations and YREC will 
search within the chosen parameter space for a solution. This 
is especially convenient since the mixing length parameter 
and in some cases the helium abundance used to compute 
stellar models must first be established from calibrated solar 
models. The calibrated values are sensitive to the choice of 
opacity tables, the equation of state formulation, the inclu- 
sion of diffusion, and the choice of model atmosphere. 



7.1 Calibrated solar models 

To produce a calibrated solar model the user inputs the age 
of the Sun and its primordial composition, i.e., mass frac- 
tion mixture of hydrogen, helium, and metals on the zero 
age main-sequence. In addition the user specifies the toler- 
ances for the luminosity and radius. YREC will then vary 
the initial value for the helium abundance and mixing length 
parameter until it has produced a model at the age of the Sun 
that has the observed radius, 6.9 58 x 10 1Q cm, and ob served 
luminosity, 3.85 15 x 10 33 erg/s dEdmonds etal.ll992h within 
the specified tolerances. When including the effects of metal 
and helium diffusion, the user has the option of inputting the 
Z/X at the surface and its allowed tolerance. In this case, 
YREC will adjust the initial helium abundance, metal abun- 
dance, and mixing length until a model at the age of the 
Sun is produced that matches the Sun's luminosity, radius, 
and surface Z/X within the specified tolerances. With 64- 
bit floating point numbers, YREC can compute a tuned solar 
model with tolerances of 1 part in 10 6 for radius and lumi- 
nosity and 1 part in 10 4 in Z/X after about 10 to 12 itera- 
tions. 

The actual procedure begins by computing one reference 
run, one run with slightly changed helium abundance, fol- 
lowed by one run with slightly changed mixing length pa- 
rameter, and then one run, if chosen, with slightly changed 
metal abundance. The luminosity, radius, and, if chosen, sur- 
face Z/X, of the final models are used to compute the deriva- 
tive matrix of luminosity, radius and surface Z/X with re- 
spect to helium abundance, mixing length parameter, and Z. 
The first order corrections to each parameter are determined 
from the derivative matrix and a new model is computed. 
The process is iterated until the model falls within the spec- 
ified tolerances. 
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7.2 Calibrated stellar models 

The process is slightly different for stars because normally 
the age of a star is unknown. Only the luminosity and surface 
temperature are used to constrain the model. In the case of 
stars, YREC adjusts the mass and either the mixing length 
parameter or the helium abundance in an attempt to produce 
a stellar evolutionary track that passes though the tolerance 
box in the theoretical HR-diagram. 

To produce a stellar model the user inputs the metal abun- 
dance, and either the helium abundance or the mixing length 
parameter. In addition the user specifies the luminosity, ef- 
fective temperature, and their corresponding tolerances. The 
user has the option of allowing the code to adjust either the 
mixing length parameter or the helium abundance. The code 
generates tracks varying the mass and the chosen parame- 
ter using a derivative matrix to produce a model that passes 
through the specified location in the HR-diagram. Once the 
optimum parameters are determined, the code computes the 
track a second time but stops the evolution when the model 
hits the specified location in the HR-diagram. The tuned 
model is constrained in mass and age. 



7.3 Pulsation models 

The pulsation output files in YREC are tai lored for the JIG 
non-adiabatic oscillation code of Guenther (lGuentherill994h . 
These files can be saved for specified models in an evolu- 
tionary sequence (say for a calibrated solar model, or for 
a calibrated stellar model), or any model for specified ages 
along the evolutionary track. 

One of the first things stellar modelers realized when us- 
ing their solar models for pulsation analysis is that the op- 
timum distribution of shells within the model for structure 
and evolutionary calculations is different from the optimum 
distribution of shells for pulsation analysis. Whereas evo- 
lutionary models need to resolve well the nuclear burning 
regions, pulsation models need to resolve the surface layers 
(for acoustic modes). For example, for evolutionary models 
great care is needed to fully resolve the thinning hydrogen 
burning shell (0.001 M© to 0.0001 M ) as the models evolve 
up the giant branch. For pulsation models it is the low den- 
sity regions, where the sound waves have the largest ampli- 
tudes, that need to be well resolved. Therefore, in order to 
produce viable models for pulsation analysis, the user in- 
creases the resolution of shells in the envelope, atmosphere, 
and the region below the base of the convection zone. Ulti- 
mately, in order to achieve frequency accuracies of the or- 
der of 1 part in 10 4 using a first order numerical pulsation 
program one needs approximately 600 shells in the interior, 
600 shells in the envelope defined as the outer 1 — 5% of the 
mass, and 600 shells in the atmosphere. To maximize self 
consistency all thermodynamic variables and their deriva- 
tives are obtained directly from the structure model. 

Related to shell resolution is the distribution of shells 
near the core. Stellar evolutionary codes do not locate a shell 



at the center owing to divide by zero complications but set 
the innermost shell a small distance away from the center. 
In order to do accurate pulsation analysis of g-modes or to 
study the p-mode small spacing parameter, both of which are 
sensitive to the structure of the deep interior, it is necessary 
to extend the innermost shell closer to the center than nor- 
mally required by stellar evolutionary calculations: compare 
1.0 x 10~ 3 radius fraction for stellar evolution to 2.0 x 10~ 7 
for stellar pulsation. A stellar model output for pulsation 
runs from the "central-most" interior shell to the top of at- 
mosphere computation near optical depth T = 10~ 10 . 



7.4 Model grids 

A useful feature of YREC is its ability to carry out exten- 
sive model calculations without user input. It is possible to 
generate in a single run, tens of thousands of evolutionary 
tracks, corresponding to tens of millions of models, covering 
a wide range of masses, compositions, mixing length param- 
eters, with each track tuned to their own particular numerical 
and ph ysical variables. This has enabled lGuenther & Browr] 
d2004l) to compute dense grids of stellar models for pulsa- 
tion analysis throughout the HR-diagram. For other grids of 
evolutionary sequences, see Sect. 1 131 



7.5 Backwards compatibility 

All new physics (e.g., opacities, nuclear reaction rates) have 
been implemented along with existing physics so that the 
user can, at any time, run YREC using older physics. 



8 Convection 

By default the local Schwarzschild criterion is implemented 
in order to determine if a mass shell is labeled as connective 
or radiative. The Ledoux stability criterion can also be used 
when a specific parameter choice is made in the local limit 
of the non-local convection treatment described below. The 
abundance of chemical species in convective cores is treated 
under the assumption of instantaneous mixing. 



9 Core Overshoot 

Since the eddy velocity at convective boundaries is non- 
zero, convective motions will penetrate into the radiative 
region. Two different forms of penetration are commonly 
distinguished: (a) inefficient penetration that does not alter 
the temperature gradient, termed "ove rmixing" here, and (b) 
subadiabatic penetration dZahn|[l99lh . where the convective 
heat transport is efficient enough to establish a nearly adia- 
batic temperature gradient. 

YREC offers a number of different options for treating 
overshoot (OS). All OS options have in common that mixing 



6 



P. Demarque, D. B. Guenther, L. H. Li, A. Mazumdar and C. W. Straka 



of chemical species in the OS region is instantaneous and all 
chemical species are homogenized within the extended zone. 
Due to the small characteristic time scale of convection in 
comparison to the thermal and nuclear timescales during the 
major burning stages this assumption holds to high accuracy. 

Among the OS options, two different approaches are dis- 
tinguished: (a) a parametric treatment where the OS extent 
is a multiple of the pressure scale height taken at the for- 
mal Schwarzschild boundary and (b) a physically motivated 
treatment where the OS extent is calculated from a non-loca l 
convection theory o riginally developed by iKuhfufil dl986l) 
and later extended bv lWuchterl & Feuchtingerld 19981) . In the 
latter case, the temperature gradient is calculated directly 
from the additional convection equation which is solved in 
addition to the canonical stellar structure equations at every 
time-step. 



9.1 Parametric Treatment 

The boundary of the OS zone is determined by adding a frac- 
tion (Xom °f th e pressure scale height to the boundary at the 
radius r$, determined by the Schwarzschild criterion: 



Table 1 Parameters of non-local convection theory 



'•new = r s + ao M Hp(r s ) 



(17) 



where the pressure scale height Hp(rs) is taken at the Schwarz- 
schild boundary. The temperature stratification in the OS 
zone is determined by the two options described above, ei- 
ther (a) the temperature gradient is not altered (overmixing) 
or (b) the temperature gradient is set to the adiabatic tem- 
perature gradient. For a fixed OCqm t he latter option produces 
larger convective cores (IZahnll 19911) . 



9.2 Non-local Convection 

As an alternative to the purely parametric treatmen t of OS, 
the on e-dimensional con vection theory deve loped bv lKuhfufil 
Jl986l) is implemented dStraka et all 120051) . In the frame- 
work of anelastic and diffusion-type approximations of the 
unknown correlation functions, KuhfuB derives one equa- 
tion for the turbulent kinetic energy from spherical averages 
of the first-order perturbed Navier-Stokes equations. The so- 
lution of this equation provides the extent of the convective 
core region and includes the effects of OS naturally, since 
the velocity of convective motions is zero where the turbu- 
lent kinetic energy vanishes. In addition, this equation also 
provides the temperature gradient in the OS region. 

9.2.1 Implemented Equations 

The new equation for the turbulent kinetic energy ft) that is 
solved in YREC is given by: 

Dft) V ad . C D _ 3/2 d , 2 



j t = -4nr 2 p 2 a t A (o 



1/2 



dm 



(18) 
(19) 



parameter 


canonical value 


description 


«MLT 


1.5 


mixing length 


a, 


0.408 


turbulent driving 


CD 


2.177 


dissipation efficiency 


a, 


0.610 a s 


overshoot distance 




1.0 


geometric mixing length 



where D/D? is the Lagrangian time derivative, p density, r 
stellar radius, m the Lagrangian mass coordinate. A, defined 
as 1/A = 1/(«mlt^p) + l/(A-f) is the geometrically lim- 
ited mixing length scale with V ac j being the adiabatic gradi- 
ent. Note that the limiting of the pressure scale height in the 
central part influences the total core size within the frame- 
work of non-local convection theories. 

A linear implicit extrapolation method is used in order to 
calculate the stationary solution of Equation (118b , i.e. Dft) /Dt : 
0. The solution yields the turbulent kinetic energy ft) at every 
mass shell. We define shells to be convective, if: 



x a < 0.1 
where 



1 



l+Fa 1 / 2 



F = 



3a s Kp 2 ACp 
I60T 3 



(20) 



(21) 



with the usual notation for the opacity fc, temperature T, spe- 
cific heat at constant pressure Cp and the Stefan-Boltzmann 
constant <7. The boundaries are sharply defined by an ex- 
treme falloff of ft) which is encountered in interior solutions 
of Equation d 1 8b - Finally, the temperature gradient can be 
calculated from: 



'ad J 



(1 



<9ft) 
an 



with 

= «t Hp f inr 2 p \ 
a s C P T\ Am J 



(22) 

(23) 
(24) 



where Am is the mass enclosed in one shell and = d In fj./d In P. 
In the case of an ideal gas with radiation pressure the di- 
mensionless parameters <5 and (p take on the values 8 = 
(4 - 3/3) //3, /3 = Pgas/P and (p = 1 respectively. In con- 
vective core regions, the temperature gradient remains very 
close to the adiabatic one whenever Xq, < 0.1. A more de- 
tailed discussion of the implemented equations and the nu- 
merica l techniques employed can be found in IStraka et alJ 
(120051) . 



9.2.2 Non-local parameters 

The implemented non-local convection theory contains five 
parameters (Table [TJ. These parameters must be calibrated, 
preferably on a well selected set of open clust ers, or on se- 
lected asteroseismic target stars like Procyon A dStraka et alJ 
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l2005h . Two of the canonical values given in Table Q] i.e. a s 
and a t can be derived by matching the KuhfuB result to mix- 
ing length theory (MLT) in the local limit (Ot = 0). 

The mixing-length parameter plays a minor role in the 
core regions, where superadiabaticity of temperature gradi- 
ent is tiny. The parameter that controls the OS zone is given 
by at and is thus the most important one to calibrate. KuhfuB 
derives Oft = 0.610 0£ s from theoretical arguments. 

In the strictly local limit (a t = 0) the KuhfuB treatment is 
equivalent to the MLT equations when based on the Ledoux 
stability criterion. 



10 Turbulence 

A method to incorporate the effects of turbulence into the 
outer layers of one-dim ensional (ID) s tellar models has been 
implemented in YREC dli et al.l2002l) . The method requires 
a detailed three-dimensional hydrodynamical simulation (3D 
RHP ) of the atmosphere an d highly superadiabatic layer of 
stars dRobinson et alj|2003h . 

The basic idea is to extract from the velocity field of the 
3D simulation three important quantities: the turbulent pres- 
sure, the turbulent kinetic energy and the anisotropy of the 
flow. Implementation into a ID stellar model thus requires 
two additional parameters, i.e. %, the specific turbulent ki- 
netic energy, and y which reflects the flow anisotropy. These 
parameters, which modify the hydrostatic equilibrium equa- 
tion and the internal energy equation, must be introduced in 
a thermodynamically self-consistent way. As a result, they 
also change the adiabatic and convective temperature gradi- 
ents, as well as the energy conservation equation. 

The next section (Sect. flOTt describes the calculation of 
X and y from the velocity field in the 3D simulation. The in- 
troduction of the parameters % and y into the stellar structure 
equations and YREC is summarized in Sect. ll0.2l Sect. 110.31 
and Sect. 110.41 The effects on /?-mode frequencies in a solar 
model are illustrated in Sect. |10.51 

10.1 Turbulent velocities 



include magnetic fields in calculating the convective temper- 
ature gr adient within the ML T framework, and used success- 
fully by iLvdon etail dl996h to explain the variation of solar 
/?-modes with the solar cycle. 

Turbulence can be measured by the turbulent Mach num- 
ber ^# = v"/v s , where v" is the turbulent velocity, and v s is 
the sound speed. The MLT is valid when ^ is sufficiently 
small. In th e outer layers of a s tar like the sun ^ can be of 
order unity dCox & Giulill 19681) . but in the deep convection 
region ^ is almost zero. The turbulent velocity is defined 
by the velocity variance: 

vr = (V-v?) 1/2 , (25) 

where the overbar denotes a combined horizontal and tem- 
poral average, and v, is the total velocity. 

Using we can define the turbulent kinetic energy per 
unit mass % as 

X = l -Jt 2 v]. (26) 

The turbulent contribution to the entropy is 

S bub =x/T, (27) 

where T is the gas temperature. 

Turbulence in the stratified layers of a stellar convection 
zone is not isotropic. We define the parameter 7 to reflect the 
anisotropy of turbulence, 

/ , tu 1 b = (y-i)pz, (28) 

where p% is the turbulent kinetic energy density. Since P tur b = 
pv" 2 , yean be related to the turbulent velocity as follows: 

y= l+2(v^/v") 2 . (29) 
7= 5/3 when turbulence is isotropic (v" = v" = v'l); 7=3 
or y = 1 when turbulence is completely anisotropic (v" = v" 
or v'l = 0, respectively). The physical meaning of 7 is the 
specific heat ratio due to turbulence. This affects the distri- 
bution of the radial turbulent pressure which is then scaled 
with the gas pressure, P gas . The total pressure is defined as 

Pt = ^gas + ^rad + ^turb ■ (30) 



The physics (thermodynamics, the equation of state, and opac- 
ities) in the 3D simulation is the same as in the ID stellar 
models. T hese simulations foll ow closely the approach de- 
scribed by lKim&Chanldl998l). and are described in more 
detail in the papers of lRobinson et"aL1 (I2003L [2004h . The full 
hydrodynamical equations were solved in a thin subsection 
of the stellar model, i.e. a 3D box located in the vicinity 
of the photosphere. For the radiative transport, the diffusion 
approximation was used in the deep region (t > 10 3 ) of the 
simulation, while the 3D Eddington approximation was used 
dUnno & St>iegellll966l) in the region above. After the simu- 
lation had reached a steady state, statistical integrations were 
performed for each simulation for over 2500 seconds in the 

case of the solar surface convection; 

For the derivation of % and 7, iLi et al. ( 2002 ) u se a s elf- 
consistent approach introduced by ILvdon & Sofial d 19951) to 



10.2 Convective temperature gradients with the turbulent 
velocities 

Since the parameters % and 7 now appear in the equation 
of state, they must be included as independent variables in 
evaluating the density derivative. We have therefore: 

dp/p = plcIPt/Pt - n'dT/T - vdx/x - v'dy/y, (31) 
where 

J" 



V = - 



dlnp \ 
d lnp \ 



As a result, the stability criterion against convection is 
modified. For similar reasons, both the convective and adia- 
batic gradients are also modified by turbulence. 
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Fig. 1 P-mode frequency difference diagrams. Turbulent model minus 
standard model (ssm), for the turbulent pressure solar model (psm), and 
the solar model with the turbulent pressure and turbulent kinetic energy 
(esm). The difference between psm and ssm is of the order of 1/iHz, 
while at high frequencies esm and ssm differ by more than 10 jiiHz. 
Plotted are the / = 0, 1, 2, 3, 4, 10, 20, . . ., 100 p-modes. 



10.3 Solar model with turbulent pressure alone 

The simplest way to take into account turbulence in solar 
modeling is to include turbulent pressure (or Reynolds stress) 
alone. In this case, only the hydrostatic equilibrium equation 
needs to be modified as follows: 
dP 



GM r 
~dM r ~ ~4nr 4 
where P = P„- ds 



(1 + J8), 



(32) 



■ rad 



P = 



2P, 



turb 



dP 



turb 



Pgr 



dP 



and 



dP 



turb 



dP 



Here 2P Va ^ / {pgr) originates from the spherical coordinate 
system adopted, representing a kind of geometric effect. The 
equations that govern the envelope integrations also need to 
be changed accordingly. 

One can construct a calibrated nonstandard model in the 
same way as one obtains the standard solar model, assuming 
that P t urb> set equal to its value for the present sun, does not 
change from the ZAMS to the present age of the sun. The 
p-mode oscillation spectrum of this calibrated solar model 
(psm) is discussed in Sect. 110.51 



10.4 Solar model with % and 7 as independent variables 

The form of the continuity equation and of the equation of 
transport of energy by radiation are not affected by turbu- 
lence. The hydrostatic equation includes a Reynolds stress 
term due to turbulence 

dP 



GM T 



1 d 

r 2 dr 



(r pv r v r ), 



(34) 



Fig. 2 P-mode frequency difference diagrams, observation minus 
model, scaled by the mode mass Q„ h for the standard solar model 
(ssm), the turbulent pressure solar model (psm), a solar model with 
fixed turbulent pressure and kinetic energy (esml), and a solar model 
with evolutionary turbulent pressure and kinetic energy (esm2, almost 
overlaps with esml). Plotted are the / = 0, 1, 2, 3, 4, 10, 20, . . ., 100 
p-modes. 



where P = P gas +/ > la d- Since the last term can be rewritten as 
dP t mb/dr + 2(Y— l)x/ r > this equation becomes 
dP T 



GM r 2(y-l)x 



(35) 



dM r 4nr 4 4nr 3 
The last term on the right hand side of Eq. (1351 also em- 
bodies the same spheric geometric effect as 2P tur b / {pgr) in 
Eq. (|33]l. 

The energy conservation equation is also modified by 
turbulence because the first law of thermodynamics must 
now include the turbulent kinetic energy. We have then: 



(33) dL r 



dM r 
where 



dt 



TdS T = c v dT - —dP r + 1 



Pru'v 



Prli'v' 



P V PHX J PHY 
The equation of energy transport by convection, 

dT T GM,- 



(36) 



rfy.(37) 



(38) 



dM r P T 4xr* conv ' 
does not change in form, but the convective temperature gra- 
dient, discussed in a previous section, is different from that 
without turbulence. The equations that govern envelope inte- 
grations also need to be changed accordingly. The oscillation 
properties of the calibrated solar model constructed under 
this assumption (esm) are discussed in the next section. 

10.5 Frequency corrections to solar p-modes 

Implementing the effects of turbulence in the outer layers of 
the stellar model modifies the calculated p-mode frequen- 
cies at high frequencies. The magnitude of the frequency 
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correction is illustrated in Fig. [Hand Fig. [2 ] for the c ase of a 
solar model, ta ken from the work of iLi et all d2002h . In the 
iLi et alJ (l2002h paper, the />mode frequencies for two cal- 
ibrated solar models that include the effects of turbulence 
are compared to the standard solar model (ssm) p-mode fre- 
quencies. The psm model is obtained by including turbulent 
pressure alone in the solar modeling, while the esm mod- 
els are obtained by introducing the turbulent variables % and 
y which include both turbulent pressure and kinetic energy. 
Contrary to a frequently made assertion, the inclusion of tur- 
bulent pressure in the pressure term has only a small effect 
on the calculated /?-mode frequencies. On the other hand, 
the inclusion of turbulent kinetic energy is significant. This 
is illustrated in Fig. Q] which shows that the frequency dif- 
ferences caused by turbulent kinetic energy are much larger 
in size than those caused by turbulent pressure alone. Fig. [2] 
indicates that the frequency changes caused by turbulent ki- 
netic energy make the computed model frequencies match 
the solar data better t han the ssm model. This result is consis- 
tent with the work of lRosenthal et al.l ( 19991) who "patched" 
a modified 3D RHD simulation by lStein & Nordlundl dl989h 
onto a ID solar model (see their Figures 1 and 6). 



10.6 SAL peak shift 

While it is always preferable to extract the y-% data from 
a 3D RHD simulation that corresponds to exactly the same 
atmospheric conditions (log g, log T e ff) as in the ID model, 
it is of interest to estimate the turbulence effects in stellar 
models where the 3D RHD simulation is not available. 

In such situations, the y-% data cannot be used directly, 
instead the data must be shif ted in order to be a pplied at the 
correct depth of the model dStraka et al.ll2006h . This shift- 
ing is motivated by an expected characteristic found in all 
3D RHD atmosphere simulations: namely that the SAL peak 
closely coincides with the turbulent pressure peak. 



Another example, in which a proper calibration is possi- 
ble with more asteroseismic data, is the detached binary sys- 
tem a Centauri: The masses of both components are known 
to high precision, the radius of the A component is mea - 
sured with interferometric techniques dKervella et al.ll2003h 
and the luminosity is also well determined through parallax 
measurements. With the help of future asteroseismic data of 
the low order />mode frequency spectrum, the stellar age of 
a Centauri can be effectively determined. Under such cir- 
cumstances the methods described to include turbulence in 
YREC are fully applicable. 



11 Seismic diagnostics 

Stellar models constructed with YREC have been used to de- 
velop seismic diagnostics to explore internal structural prop- 
erti es of stars that cou ld not be observed by any other means. 

Bas u et al. I d200l have used low degree acoustic modes 
to determine the helium abundance in the envelopes of low- 
mass main sequence stars with precision. The oscillatory 
signal in the frequencies caused by the depression in Fi in 
the second helium ionization zone is used. For frequency er- 
rors of one part in 10 4 , the expected ay in the estimated Y 
ranges from 0.03 for O.8M stars to 0.01 for 1 .2M Q main se- 
quence stars. In more evolved stars, this approach is feasible 
if t he relative errors i n the f requencies are less than 10~ 4 . 

iMazumdar et all (120061) have explored asteroseismic di- 
agnostics of convective core mass using small frequency sep- 
arations of low-degree p-modes. Small separations can also 
be combined to derive convective core overshoot diagnos- 
tics. It was shown that in stars with convective cores, the 
mass of the convective core can be estimated to within 5% 
if the total mass is known, although systematic errors in the 
total mass could introduce errors of up to 20%. The evolu- 
tionary stage of the star, determined in terms of the central 
hydrogen content is much less sensitive to the mass estimate. 



10.7 Calibration 

The presented method for including turbulence does not re- 
move MLT and therefore the uncertainties inherent in the 
mixing length parameter remain. In order to make quantita- 
tive predictions, both mass and age of the star must be known 
to high precision in addition to the luminosity L and effec- 
tive temperature T e ^. Instead of the latter an interferometric 
radius measurement is preferable, since the measurement is 
usually more precise. 

In the case of the Sun, the age is known to high precision 
and the mixing length parameter and hydrogen mass frac- 
tion can be calibrated to the known solar lumi n osity and ra- 
dius. A s demonstrated in lGuenther et al.l d2005h : ls"traka et all 
d2006h . asteroseismic data can be instrumental in other stars, 
since low order /?-mode frequencies anchor the interior model 
effectively in age and mass. When calibrated to the same lu- 
minosity and effective temperature, a differential assessment 
of turbulence effects can be derived. 



12 Solar neutrinos and helioseismology 

Different versions of YREC have been used to study the 
structure of the sola r interior, the solar neutrino problem 
and helioseismology dGuenther & Demarqudl 19971) and ref- 
er ences therein. In part i cular, the important series of papers 
bv lBahcall et"ai1 (12004 l2005h on solar neutrinos, helioseis- 
mology and solar abundances also made use of a dedicated 
version of YREC. 



13 Other YREC applications 

A variety of applications to stellar structure theory and evo- 
lution have been carried out using YREC. In addition of the 
wor k on stellar rotat i on me ntio ned in the intro duction (see 
also IChaboyer etail dl995h and iBarnesI d2003h l. one notes 
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the pi oneer work on stellar collisions and mergers dSills et al.1 
1 19971) . 

Important research in stellar population studies and pop- 
ulation synthesis continues to be c arried out with the Yonsei- 
Yale isochrones (YY isochrones) dYi et al.l200lh . Frequently 
quoted research on helium burning phases of evoluti on (hori- 
zontal-branch) has a lso been carried out with YREC dLee et all 
ll994l : lYietalll997h . 
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